#The following script generates the changepoint analysis and corresponding plots (Figures 2 and A6 in the article and Supplemental Appendix)
#It uses the following files: weekly_cyclone_ts.dta and changepoint-util.R

#read-in or source-in necessary packages
library(MCMCpack)
library(MASS)
library(foreign)
source("changepoint-util.R")

#read-in weekly data
cyclonedata <- read.dta(file = "weekly_cyclone_ts.dta")

#identify the number of time points for the data
num.days <- nrow(cyclonedata)

#######################################
####Main Negative Binomial Approach####
#######################################

#setup model parameters
mcmc <- 100000
burnin <- 5000
thin <- 100
v <- 10000
m <- 15

## create smart starting values for rho/P
hh <- cut(1:nrow(cyclonedata), seq(1,nrow(cyclonedata), 40))
rho.start <- rep(summary(glm.nb(moz_rebcit_acled ~ 1, data = cyclonedata))$theta, times = m)
P.start <- matrix(0.05/(m-1), nrow = m, ncol = m)
diag(P.start) <- rep(0.95, times = m)

out <- HDPHMMnegbin(moz_rebcit_acled ~ 1, data =cyclonedata,
                    K = m, verbose = v,
                    mcmc = mcmc, burnin = burnin, thin = thin,
                    a.alpha = 1, b.alpha = 0.1,
                    a.gamma = 1, b.gamma = 0.1,
                    a.theta = 100, b.theta = 1,
                    e = 2, f = 2, g = 10,
                    b0=0, B0 = 1/9, P.start = P.start,
                    rho.step = 0.75,
                    rho.start = rho.start, seed = list(NA, 2),
                    theta.start = 0.95, gamma.start = 10, ak.start = 10)

post.mean <- post.rho <- rep(NA, times = nrow(cyclonedata))
post.rho <- matrix(NA, nrow = nrow(out), ncol = nrow(cyclonedata))
post.beta <- matrix(NA, nrow = nrow(out), ncol = nrow(cyclonedata))
for (j in 1:nrow(cyclonedata)) {
  inds <- c(1:nrow(out)) + (nrow(out)) * (attr(out, "s.store")[,j]-1)
  post.beta[,j] <- out[inds]
  post.rho[,j] <- attr(out, "rho.store")[inds]
}
ch.pt <- rbind(0, diff(t(attr(out, "s.store"))) != 0)
pr.st <- rowMeans(ch.pt)
post.mean <- colMeans(post.beta)
nb <- names(which.max(table(attr(out, "num.regimes"))))
nb <- as.numeric(nb) - 1
ints <- find.intervals(pr.st, prob = 0.95)

fulldata <- cyclonedata[, c("year","week", "moz_rebcit_acled",
                          "fullcyclone")]
full.data <- cbind(fulldata, pr.st, post.mean)

positivetrans<-0
for(i in 1:nrow(full.data)){
positivetrans<-c(positivetrans,(full.data$moz_rebcit_acled[i]-full.data$moz_rebcit_acled[i-1]))
}
full.data$positivetrans<-positivetrans
full.data$pr.st.post<-ifelse(full.data$positivetrans>0,full.data$pr.st,0)
full.data$pr.st.neg<-ifelse(full.data$positivetrans<0,full.data$pr.st,0)
full.data$pr.st.large<-ifelse(full.data$pr.st>0.5,full.data$pr.st,0)

cumcounts<-0
for(i in 1:nrow(full.data)){
cumcounts<-c(cumcounts,(cumcounts[i]+full.data$moz_rebcit_acled[i]))
}

#Plot Figure 2
par(mar = c(5.1, 4.1, 2.1, 4.1))
plot(x = as.numeric(rownames(full.data)), y = full.data$pr.st.large, type = "h", lwd = 2,
     col = rgb(0,0,0, alpha = 0.5),
     ylim = c(0,1), bty = "n", yaxt = "n",
     xaxt = "n", ylab = "Probability of a Changepoint",xlab="Weekly Observations: March 2018 - April 2020")
abline(v=55,col="red",lwd=3)
abline(v=61,col="blue",lwd=3)
text(49.5,.9,"Cyclone Idai",cex=.85,col="red")
text(68.5,.9,"Cyclone Kenneth",cex=.85,col="blue")
axis(1, at=c(0, 20, 40, 60, 80, 100),labels=c("02/2018", "07/2018", "12/2018", "04/2019", "09/2019", "01/2020")) 
axis(2, at=c(0.0, 0.2, 0.4, 0.6, 0.8, 1.0),labels=c("0.0","0.2", "0.4", "0.6", "0.8", "1.0")) 
box()
par(new=TRUE)
plot(x = as.numeric(rownames(full.data)),full.data$moz_rebcit_acled,xlab="",main="",ylab="",yaxt="n",xaxt="n")
axis(4, at=c(0,2,4,6,8,10),labels=c("0","2", "4", "6", "8", "10")) 
mtext("Insurgent Violence Against Civilians", side = 4, line =2.7)


###################################
####Poisson Robustness Approach####
###################################

mm <- apply(attr(out, "s.store"), 1,
            function(x) sum(diff(x) != 0, na.rm = TRUE))
mmm <- median(mm)

out.po <- MCMCpoissonChange(moz_rebcit_acled ~ 1, data =cyclonedata, m = 13,
                            marginal.likelihood = "Chib95", verbose = v,
                            mcmc = mcmc, burnin = burnin, thin = thin,
                            a = 20, b = 0.1,
                            ##e = 2, f = 2, d = 1,
                            b0 = 0, B0 = 1/100, c0=2, d0=1)
po.beta <- matrix(NA, nrow = nrow(out.po), ncol = nrow(cyclonedata))
for (j in 1:nrow(cyclonedata)) {
  inds <- c(1:nrow(out.po)) + (nrow(out.po)) * (attr(out.po, "s.store")[,j]-1)
  po.beta[,j] <- out.po[inds]
}
ch.pt <- rbind(0, diff(t(attr(out.po, "s.store"))) != 0)
pr.st.po <- rowMeans(ch.pt)
po.mean <- colMeans(po.beta)
full.data$pr.st.po <- pr.st.po
full.data$po.mean <- po.mean
full.data$pr.st.po<-ifelse(full.data$pr.st>0.5,full.data$pr.st.po,0)

#Plot Figure A.6
par(mar = c(5.1, 4.1, 2.1, 4.1))
plot(x = as.numeric(rownames(full.data)), y = full.data$pr.st.po, type = "h", lwd = 2,
     col = rgb(0,0,0, alpha = 0.5),
     ylim = c(0,1), bty = "n", yaxt = "n",
     xaxt = "n", ylab = "Probability of a Changepoint",xlab="Weekly Observations: March 2018 - April 2020")
abline(v=55,col="red",lwd=3)
abline(v=61,col="blue",lwd=3)
text(49.5,.9,"Cyclone Idai",cex=.85,col="red")
text(68.5,.9,"Cyclone Kenneth",cex=.85,col="blue")
axis(1, at=c(0, 20, 40, 60, 80, 100),labels=c("02/2018", "07/2018", "12/2018", "04/2019", "09/2019", "01/2020")) 
axis(2, at=c(0.0, 0.2, 0.4, 0.6, 0.8, 1.0),labels=c("0.0","0.2", "0.4", "0.6", "0.8", "1.0")) 
box()
par(new=TRUE)
plot(x = as.numeric(rownames(full.data)),full.data$moz_rebcit_acled,xlab="",main="",ylab="",yaxt="n",xaxt="n")
axis(4, at=c(0,2,4,6,8,10),labels=c("0","2", "4", "6", "8", "10")) 
mtext("Insurgent Violence Against Civilians", side = 4, line =2.7)
